
meff.m <- function(form, allmods, dums, newdat, reps=1000){
#pull data
    nk <- length(dums)
    vnames <- all.vars(form[[3]])
    if (vnames[2]=="j") vnames<- vnames[-2]
    vtank  <- matrix(NA, m*n, k-1) #holding matrix for variables
    
    for (ia in 1:(k-1)){
        for (ib in 1:m){
            if(ia>1){ 
                datnow <- newdat$imputations[[ib]]
                attach(datnow)
                vtank[(1+(ib-1)*n):(n+(ib-1)*n), ia] <- get(paste(vnames[ia]))
                detach(datnow)
            } else 
            vtank[(1+(ib-1)*n):(n+(ib-1)*n), ia] <- imr[,ib]
        } 
    }
    xprof <- c(1,apply(vtank, 2, median))
    meff.aux <- array(NA, c(3, k,m))
    meff90.aux <- array(NA, c(3, k,m))
    meff.select.full <- matrix(0, 3, k)
    meff90.select.full <- matrix(0, 3, k)
    coef.select.full <- rep(0, k)
    for(j in 1:m){
        aux.mod<-allmods[[j]]
        varcov <- sandwich(aux.mod)
        scoef <- t(mvrnorm(reps, aux.mod$coef, varcov))   #draws from sampling distribution
        coef.select.full <- coef.select.full+aux.mod$coef
        for(i in 1:nk){
            xaux2 <- xprof
            if (dums[i]==1 & i>1) xaux2[i] <- 0 #set dummies to 0
            marg.raw <- dnorm(t(as.matrix(xaux2))%*%scoef)*scoef[i,]
            meff.aux[,i,j] <- quantile(marg.raw, probs=c(.025, .5, .975))
            meff90.aux[,i,j] <- quantile(marg.raw, probs=c(.05, .5, .95))
        }
            meff.select.full <- meff.select.full + meff.aux[,,j]
            meff90.select.full <- meff90.select.full + meff90.aux[,,j]
    }
    meff.select.full <- meff.select.full/newdat$m
    meff90.select.full <- meff90.select.full/newdat$m
    coef.select.full<- coef.select.full/newdat$m
    #on the basis of a predicted probability of agreement of 
    eofx.select.full <- pnorm(t(as.matrix(xprof))%*% coef.select.full)
    colnames(meff.select.full) <- c("const", vnames)
    colnames(meff90.select.full) <- c("const", vnames)
    list(meff=meff.select.full, meff90=meff90.select.full, coef=coef.select.full, exb=eofx.select.full)

}

meff.m.int <- function(form, allmods, dums, newdat, reps=1000, pos){
#pull data
    nk <- length(dums)
    vnames <- all.vars(form[[3]])
    if (vnames[2]=="j") vnames<- vnames[-2]
    vtank  <- matrix(NA, m*n, k-1) #holding matrix for variables
    
    for (ia in 1:(k-1)){
        for (ib in 1:m){
            if(ia>1){ 
                datnow <- newdat$imputations[[ib]]
                attach(datnow)
                vtank[(1+(ib-1)*n):(n+(ib-1)*n), ia] <- get(paste(vnames[ia]))
                detach(datnow)
            } else 
            vtank[(1+(ib-1)*n):(n+(ib-1)*n), ia] <- imr[,ib]
        } 
    }
    xprof <- c(1,apply(vtank, 2, median))
    xprof[pos] <- 0
    meff.aux <- array(NA, c(3, k,m))
    meff90.aux <- array(NA, c(3, k,m))
    meff.select.full <- matrix(0, 3, k)
    meff90.select.full <- matrix(0, 3, k)
    coef.select.full <- rep(0, k)
    for(j in 1:m){
        aux.mod<-allmods[[j]]
        varcov <- sandwich(aux.mod)
        scoef <- t(mvrnorm(reps, aux.mod$coef, varcov))   #draws from sampling distribution
        coef.select.full <- coef.select.full+aux.mod$coef
        for(i in 1:nk){
            xaux2 <- xprof
            if (dums[i]==1 & i>1) xaux2[j] <- 0 #set dummies to 0
            marg.raw <- dnorm(t(as.matrix(xaux2))%*%scoef)*scoef[i,]
            meff.aux[,i,j] <- quantile(marg.raw, probs=c(.025, .5, .975))
            meff90.aux[,i,j] <- quantile(marg.raw, probs=c(.05, .5, .95))
        }
            meff.select.full <- meff.select.full + meff.aux[,,j]
            meff90.select.full <- meff90.select.full + meff90.aux[,,j]
    }
    meff.select.full <- meff.select.full/newdat$m
    meff90.select.full <- meff90.select.full/newdat$m
    coef.select.full<- coef.select.full/newdat$m
    #on the basis of a predicted probability of agreement of 
    eofx.select.full <- pnorm(t(as.matrix(xprof))%*% coef.select.full)
    colnames(meff.select.full) <- c("const", vnames)
    colnames(meff90.select.full) <- c("const", vnames)
    meff.select.full.0 <- meff.select.full
    meff90.select.full.0 <- meff90.select.full
    
    xprof[pos] <- 1
    meff.aux <- array(NA, c(3, k,m))
    meff90.aux <- array(NA, c(3, k,m))
    meff.select.full <- matrix(0, 3, k)
    meff90.select.full <- matrix(0, 3, k)
    coef.select.full <- rep(0, k)
    for(j in 1:m){
        aux.mod<-allmods[[j]]
        varcov <- sandwich(aux.mod)
        scoef <- t(mvrnorm(reps, aux.mod$coef, varcov))   #draws from sampling distribution
        coef.select.full <- coef.select.full+aux.mod$coef
        for(i in 1:nk){            
            xaux2 <- xprof
            if (dums[i]==1 & i>1) xaux2[j] <- 0 #set dummies to 0
            marg.raw <- dnorm(t(as.matrix(xaux2))%*%scoef)*scoef[i,]
            meff.aux[,i,j] <- quantile(marg.raw, probs=c(.025, .5, .975))
            meff90.aux[,i,j] <- quantile(marg.raw, probs=c(.05, .5, .95))
        }
            meff.select.full <- meff.select.full + meff.aux[,,j]
            meff90.select.full <- meff90.select.full + meff90.aux[,,j]
    }
    meff.select.full <- meff.select.full/newdat$m
    meff90.select.full <- meff90.select.full/newdat$m
    coef.select.full<- coef.select.full/newdat$m
    #on the basis of a predicted probability of agreement of 
    eofx.select.full <- pnorm(t(as.matrix(xprof))%*% coef.select.full)
    colnames(meff.select.full) <- c("const", vnames)
    colnames(meff90.select.full) <- c("const", vnames)
    meff.select.full.1 <- meff.select.full
    meff90.select.full.1 <- meff90.select.full
    
    
    list(meff.0=meff.select.full.0, meff90.0=meff90.select.full.0, 
            meff.1=meff.select.full.1, meff90.1=meff90.select.full.1, 
            coef=coef.select.full, exb=eofx.select.full)
}


meff.s<- function(mod, form, dums, newdat, reps=1000, npol){
    #npol is the number of auxiliary vars that won't have marginal effects computed
    xaux <- setx(mod, fn=list(numeric=median, ordered=median, other=mode))
    meff.aux <- array(NA, c(3, length(dums)-npol,m))
    meff90.aux <- array(NA, c(3, length(dums)-npol,m))
    meff.select.full <- matrix(0, 3, length(dums)-npol)
    meff90.select.full <- matrix(0, 3, length(dums)-npol)
    coef.select.full <- rep(0, length(dums))
    for(i in 1:m){
        aux.mod<-glm(form, 
                data=newdat$imputations[[i]], binomial(link = "probit"), na.action=na.exclude)
        scoef <- t(mvrnorm(reps, aux.mod$coef, summary(aux.mod)$cov.unscaled))   #draws from sampling distribution
        coef.select.full <- coef.select.full+aux.mod$coef
        for(j in 1:(length(dums)-npol)){
            xprof <- xaux        
            if (dums[j]) xprof[j] <- 0 #set dummies to 0
            marg.raw <- dnorm(as.matrix(xprof)%*%scoef)*scoef[j,]
            meff.aux[,j,i] <- quantile(marg.raw, probs=c(.025, .5, .975))
            meff90.aux[,j,i] <- quantile(marg.raw, probs=c(.05, .5, .95))
        }
            meff.select.full <- meff.select.full + meff.aux[,,i]
            meff90.select.full <- meff90.select.full + meff90.aux[,,i]
    }
    meff.select.full <- meff.select.full/newdat$m
    meff90.select.full <- meff90.select.full/newdat$m
    coef.select.full<- coef.select.full/newdat$m
    #on the basis of a predicted probability of agreement of 
    eofx.select.full <- pnorm(as.matrix(xaux) %*% coef.select.full)
    list(meff=meff.select.full, meff90=meff90.select.full, coef=coef.select.full, exb=eofx.select.full)
}

meff.s.int<- function(mod, form, dums, newdat, reps=1000, npol, pos){
    #npol is the number of auxiliary vars that won't have marginal effects computed
    xaux <- setx(mod, fn=list(numeric=median, ordered=median, other=mode))
    meff.aux <- array(NA, c(3, length(dums)-npol,m))
    meff90.aux <- array(NA, c(3, length(dums)-npol,m))
    meff.select.full <- matrix(0, 3, length(dums)-npol)
    meff90.select.full <- matrix(0, 3, length(dums)-npol)
    coef.select.full <- rep(0, length(dums))
    for(i in 1:m){
        aux.mod<-glm(form, 
                data=newdat$imputations[[i]], binomial(link = "probit"), na.action=na.exclude)
        scoef <- t(mvrnorm(reps, aux.mod$coef, summary(aux.mod)$cov.unscaled))   #draws from sampling distribution
        coef.select.full <- coef.select.full+aux.mod$coef
        for(j in 1:(length(dums)-npol)){
            xprof <- xaux
            xprof[pos] <- 0
            if (dums[j]) xprof[j] <- 0 #set dummies to 0
            marg.raw <- dnorm(as.matrix(xprof)%*%scoef)*scoef[j,]
            meff.aux[,j,i] <- quantile(marg.raw, probs=c(.025, .5, .975))
            meff90.aux[,j,i] <- quantile(marg.raw, probs=c(.05, .5, .95))
        }
            meff.select.full <- meff.select.full + meff.aux[,,i]
            meff90.select.full <- meff90.select.full + meff90.aux[,,i]
    }
    meff.select.full <- meff.select.full/newdat$m
    meff90.select.full <- meff90.select.full/newdat$m
    meff.select.full.0 <- meff.select.full
    meff90.select.full.0 <- meff90.select.full
 
 for(i in 1:m){
        aux.mod<-glm(form, 
                data=newdat$imputations[[i]], binomial(link = "probit"), na.action=na.exclude)
        scoef <- t(mvrnorm(reps, aux.mod$coef, summary(aux.mod)$cov.unscaled))   #draws from sampling distribution
        coef.select.full <- coef.select.full+aux.mod$coef
        for(j in 1:(length(dums)-npol)){
            xprof <- xaux
            xprof[pos] <- 1
            if (dums[j]) xprof[j] <- 0 #set dummies to 0
            marg.raw <- dnorm(as.matrix(xprof)%*%scoef)*scoef[j,]
            meff.aux[,j,i] <- quantile(marg.raw, probs=c(.025, .5, .975))
            meff90.aux[,j,i] <- quantile(marg.raw, probs=c(.05, .5, .95))
        }
            meff.select.full <- meff.select.full + meff.aux[,,i]
            meff90.select.full <- meff90.select.full + meff90.aux[,,i]
    }
    meff.select.full <- meff.select.full/newdat$m
    meff90.select.full <- meff90.select.full/newdat$m
    meff.select.full.1 <- meff.select.full
    meff90.select.full.1 <- meff90.select.full
    
    coef.select.full<- coef.select.full/newdat$m
    #on the basis of a predicted probability of agreement of 
    eofx.select.full <- pnorm(as.matrix(xaux) %*% coef.select.full)
    list(meff.0=meff.select.full.0, meff90.0=meff90.select.full.0,
        meff.1=meff.select.full.1, meff90.1=meff90.select.full.1, 
            coef=coef.select.full, exb=eofx.select.full)
}

#effect of regime
eff.s.reg<- function(mod, form, dums, newdat, reps=1000, npol, pos, number){
    #npol is the number of auxiliary vars that won't have marginal effects computed
    xaux <- setx(mod, fn=list(numeric=median, ordered=median, other=mode))
    eff.aux <- array(NA, c(3, length(dums)-npol,m))
    eff90.aux <- array(NA, c(3, length(dums)-npol,m))
    eff.select.full <- matrix(0, 3, length(dums)-npol)
    eff90.select.full <- matrix(0, 3, length(dums)-npol)
    coef.select.full <- rep(0, length(dums))
    for(i in 1:m){
        aux.mod<-glm(form, 
                data=newdat$imputations[[i]], binomial(link = "probit"), na.action=na.exclude)
        scoef <- t(mvrnorm(reps, aux.mod$coef, summary(aux.mod)$cov.unscaled))   #draws from sampling distribution
        coef.select.full <- coef.select.full+aux.mod$coef
        for(j in 1:(length(dums)-npol)){
            xprof <- xaux
            xprof[pos] <- 0
            if (dums[j]) xprof[j] <- 0 #set dummies to 0
            eff.raw <- pnorm(as.matrix(xprof)%*%scoef)
            eff.aux[,j,i] <- quantile(eff.raw, probs=c(.025, .5, .975))
            eff90.aux[,j,i] <- quantile(eff.raw, probs=c(.05, .5, .95))
        }
            eff.select.full <- eff.select.full + eff.aux[,,i]
            eff90.select.full <- eff90.select.full + eff90.aux[,,i]
    }
    eff.select.full <- eff.select.full/newdat$m
    eff90.select.full <- eff90.select.full/newdat$m
    eff.select.full.0 <- eff.select.full
    eff90.select.full.0 <- eff90.select.full
 
 for(i in 1:m){
        aux.mod<-glm(form, 
                data=newdat$imputations[[i]], binomial(link = "probit"), na.action=na.exclude)
        scoef <- t(mvrnorm(reps, aux.mod$coef, summary(aux.mod)$cov.unscaled))   #draws from sampling distribution
        coef.select.full <- coef.select.full+aux.mod$coef
        for(j in 1:(length(dums)-npol)){
            xprof <- xaux
            xprof[pos] <- number
            if (dums[j]) xprof[j] <- 0 #set dummies to 0
            eff.raw <- pnorm(as.matrix(xprof)%*%scoef)
            eff.aux[,j,i] <- quantile(eff.raw, probs=c(.025, .5, .975))
            eff90.aux[,j,i] <- quantile(eff.raw, probs=c(.05, .5, .95))
        }
            eff.select.full <- eff.select.full + eff.aux[,,i]
            eff90.select.full <- eff90.select.full + eff90.aux[,,i]
    }
    eff.select.full <- eff.select.full/newdat$m
    eff90.select.full <- eff90.select.full/newdat$m
    eff.select.full.1 <- eff.select.full
    eff90.select.full.1 <- eff90.select.full
    
    coef.select.full<- coef.select.full/newdat$m
    #on the basis of a predicted probability of agreement of 
    eofx.select.full <- pnorm(as.matrix(xaux) %*% coef.select.full)
    list(eff.0=eff.select.full.0, eff90.0=eff90.select.full.0,
        eff.1=eff.select.full.1, eff90.1=eff90.select.full.1, 
            coef=coef.select.full, exb=eofx.select.full)
}








meff.s.noUS <- function(form, dums, newdat, reps=1000){
    getmedian <- function(j, dat){
        apply(dat$imputations[[j]][dat$imputations[[j]]$usa==0,], 2, median)
    }
    rmed <- sapply(1:newdat$m, getmedian, dat=newdat)
    xaux <- apply(rmed, 1, median)
    xaux <- data.frame(t(as.matrix(xaux)))
    xaux <- model.matrix(form, xaux)
    meff.aux <- array(NA, c(3, length(dums)-12,newdat$m))
    meff90.aux <- array(NA, c(3, length(dums)-12,newdat$m))
    meff.select.full <- matrix(0, 3, length(dums)-12)
    meff90.select.full <- matrix(0, 3, length(dums)-12)
    coef.select.full <- rep(0, length(dums))
    for(i in 1:m){
        aux.mod<-glm(form, 
                data=newdat$imputations[[i]][newdat$imputations[[i]]$usa==0,], binomial(link = "probit"), na.action=na.exclude)
        scoef <- t(mvrnorm(reps, aux.mod$coef, summary(aux.mod)$cov.unscaled))   #draws from sampling distribution
        coef.select.full <- coef.select.full+aux.mod$coef
        for(j in 1:(length(dums)-12)){
            xprof <- xaux        
            if (dums[j]) xprof[j] <- 0 #set dummies to 0
            marg.raw <- dnorm(as.matrix(xprof)%*%scoef)*scoef[j,]
            meff.aux[,j,i] <- quantile(marg.raw, probs=c(.025, .5, .975))
            meff90.aux[,j,i] <- quantile(marg.raw, probs=c(.05, .5, .95))
        }
            meff.select.full <- meff.select.full + meff.aux[,,i]
            meff90.select.full <- meff90.select.full + meff90.aux[,,i]
    }
    meff.select.full <- meff.select.full/newdat$m
    meff90.select.full <- meff90.select.full/newdat$m
    coef.select.full<- coef.select.full/newdat$m
    #on the basis of a predicted probability of agreement of 
    eofx.select.full <- pnorm(as.matrix(xaux) %*% coef.select.full)
    colnames(meff.select.full) <- colnames(xaux)[1:(length(dums)-12)]
    colnames(meff90.select.full) <- colnames(xaux)[1:(length(dums)-12)]
    list(meff=meff.select.full, meff90=meff90.select.full, coef=coef.select.full, exb=eofx.select.full)
}



meff.s.dem <- function(form, dums, newdat, reps=1000){
    getmedian <- function(j, dat){
        apply(dat$imputations[[j]][dat$imputations[[j]]$jointdem==1,], 2, median, na.rm=TRUE)
    }
    rmed <- sapply(1:newdat$m, getmedian, dat=newdat)
    xaux <- apply(rmed, 1, median)
    xaux <- data.frame(t(as.matrix(xaux)))
    xaux <- model.matrix(form, xaux)
    meff.aux <- array(NA, c(3, length(dums)-12,newdat$m))
    meff90.aux <- array(NA, c(3, length(dums)-12,newdat$m))
    meff.select.full <- matrix(0, 3, length(dums)-12)
    meff90.select.full <- matrix(0, 3, length(dums)-12)
    coef.select.full <- rep(0, length(dums))
    for(i in 1:m){
        aux.mod<-glm(form, 
                data=newdat$imputations[[i]][newdat$imputations[[i]]$jointdem==1,], binomial(link = "probit"), na.action=na.exclude)
        scoef <- t(mvrnorm(reps, aux.mod$coef, summary(aux.mod)$cov.unscaled))   #draws from sampling distribution
        coef.select.full <- coef.select.full+aux.mod$coef
        for(j in 1:(length(dums)-12)){
            xprof <- xaux        
            if (dums[j]) xprof[j] <- 0 #set dummies to 0
            marg.raw <- dnorm(as.matrix(xprof)%*%scoef)*scoef[j,]
            meff.aux[,j,i] <- quantile(marg.raw, probs=c(.025, .5, .975))
            meff90.aux[,j,i] <- quantile(marg.raw, probs=c(.05, .5, .95))
        }
            meff.select.full <- meff.select.full + meff.aux[,,i]
            meff90.select.full <- meff90.select.full + meff90.aux[,,i]
    }
    meff.select.full <- meff.select.full/newdat$m
    meff90.select.full <- meff90.select.full/newdat$m
    coef.select.full<- coef.select.full/newdat$m
    #on the basis of a predicted probability of agreement of 
    eofx.select.full <- pnorm(as.matrix(xaux) %*% coef.select.full)
    colnames(meff.select.full) <- colnames(xaux)[1:(length(dums)-12)]
    colnames(meff90.select.full) <- colnames(xaux)[1:(length(dums)-12)]
    list(meff=meff.select.full, meff90=meff90.select.full, coef=coef.select.full, exb=eofx.select.full)
}
